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ABSTRACT 


The theory and background of the algebraic substitution 
synthesis method for the digital filters from continuous 
filter characterizations is presented with emphasis on the 
frequency distortion phenomenon. An analysis of the Forward 
Euler, Backward Euler and Trapezoidal numerical integration 
algorithms is undertaken and appropriate transformations 
are obtained. A general integration formula, encompassing 
the above algorithms as special cases, is analyzed and its 
application to the synthesis problem is pointed out. Direct 
transformations for discrete filter frequency response 
characteristics from continuous filter characterizations are 


derived. 
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foo LN ERODUCETRO 


Digital filter design can be accomplished using several 
different techniques. [1] One particularly convenient method 
as to transform a continuous filter sinto a digitale couse 
means of some substitution procedure. This offers the designer 
a large body of filter designs from which to draw the form 
applicable to a particular task. However, this procedure has 
certain drawbacks. Chief among these is the distortion which 
results when transforming frequencies from the continuous 
domain into frequencies in the discrete domain. An under- 
Standing of this frequency distortion is essential when the 
substitution method is employed. Fig. 1.1 shows the relation- 
ship between the continuous filter and the digital filter. 

The transfer function of the continuous filter is derived 
from a differential equation by means of the Laplace Transform. 
The frequency response of the continuous filter can be eval- 
uated by letting s = jw as shown in the upper right of the 
diagram. A difference equation can be formed by applying a 
numerical integration algorithm to the continuous differential 
equation. This difference equation can then be z-transformed 
to arrive at a discrete transfer function as shown. This 
discrete transfer function can also be developed directly by 
meplying a substitution (s= £(z)) to the continuous transter 
function. This substitution is derived from the numerical 


integration algorithm and the two discrete transfer functions 
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thus obtained are identical. This procedure is the key to the 
algebraic substitution technique. 

The frequency response of the discrete transfer function 
can now be evaluated by letting z = exp (jM). The spectrum 
thus obtained will not, in general, be identical to that of 
the continuous filter spectrum. The difference between the two 
gives rise to frequency distortion. The purpose of this thesis 
is to analyze this frequency distortion for several numerical 
integration algorithms. 

The source of the frequency distortion can be seen in 
Fig. 1.2. This is a representation of the effect of the s-plane 
to z-plane transformation on the sinusoidal steady-state fre- 
quency domain characteristics of the filter. If the transfor- 
mation was linear as depicted by the dotted line, the continuous 
spectrum would be perfectly reproduced in the discrete fre- 
quency domain. However, the transformation is not linear. Its 
shape, of course, depends on the particular integration algorithm 
used and the transforming function is periodic in nature. The 
solid line represents an arbitrary transforming function and 
its effect on the frequency characteristics of the discrete 
filter. The distortion that results from this nonlinear trans- 
formation is the subject of the analysis that follows. 

Much use is made of the Z - Transform technique to accomplish 
this analysis and the theory and use of this technique is 
described by many authors including E. I. Jury [2] and B. C. 

Kuo [3]. A paper comparing various digital integrators was 
written by J. R. Salzer [4] and is one of the classic references 


on this subject. Gold and Rader [5] have discussed the 
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substitution technique in a limited form, and Golden and 
Kaiser [6] have discussed the bilinear substitution, which 
ee to trapezoidal integration, in a thorough manner. 
Chapter II of this thesis contains the theory and develop- 
ment of the algebraic substitution technique and describes how 
a filter design can be accomplished using this method. Chapter 
III presents a detailed analysis of three specific numerical 
integration algorithms ~- Backward and Forward Euler and the 
Trapezoidal Rule. The effect of frequency distortion is shown 
and analyzed. Also specific limitations on the use of the 
Buler methods are pointed out. Chapter IV discusses a general 
integration formula (which includes the above as special cases) 
and how it can be employed. General case Sinusoidal steady- 
state frequency transformations are then derived, and an 
example of a filter synthesis is presented to demonstrate the 
distortion problem. Chapter V contains the summary and 


Suggestions for further research. 
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Tl. THEORY AND DEVELCR MEE 
OF ALGEBRAIC SUBSTITULION DESC ME TneE 

There are many instances when a discrete time representa- 
tion of a continuous system or filter is required. Digital 
Simulation of a control system for computer analysis is one 
example. Another important example is the extension of con- 
tinuous filter theory to the design and synthesis of digital 
filters. This application in particular provides incentive 
for developing a method of design which will ingoepernee the 
great wealth of information that has been accumulated. for 
continuous filters. Designs for continuous filters are well- 
known and have been tabulated in various handbooks [7]. Thus 
a direct method for obtaining a digital realization of a 
continuous filter is not only theoretically appealing but has 
definite practical benefits. 

The usual procedure for implementing a digital filter is 
to first express the filter as a transfer function in the z- 
domain. This transfer function iS normally expressed as a 
ratio of two polynomials in inverse powers of z, where zt 
represents a time delay [8]. This ratio then indicates the 
linear operations which must be performed upon the past values 
ee the input and output samples to obtain the present value 
of the output. For example, consider the following transfer 


manction 





H(Z) = (28) 


idk 





Since this transfer function is the ratio SGieene transite al 
of the output to the transform of the input this can be 


expressed as 





= (2-2) 


The latter expression is the recursive form of the digital 
filter transfer function. 

In a like manner, the most common description of continuous 
filters in the Laplace or s-domain is as the ratio of two 
polynomials ins. Consequently any substitution that can be 
used to directly replace the variable s in the continuous 
transfer function with a function of z enables a discrete 


filter realization to be formed directly. Thus if 
s = £(z) (2-3) 
Then H(s) becomes 


G(z) = H(s) | = H[f(z)] (2-4) 
s = £(z) 
Where H(s) is the continuous filter transfer function and 
f(z) is a generic algebraic expression. 
If the function f(z) is a ratio of two polynomials of 
inverse powers of z, and if H(s) is a ratio of two polynomials 
in s, then it is a simple matter to show that the G(z) thus 


obtained will always be the ratio of two polynomials of 


2 





inverse powers of z. The process expressed by Eq. (2-4) is an 
algebraic substitution method which converts a rational 
polynomial in s to a rational polynomial in z. It is the 


nature of this direct substitution which is discussed in detail. 


A. DEVELOPMENT OF THE ALGEBRAIC SUBSTITUTION 

The question that must now be answered is what substitu- 
tions can be used in Eq. (2-4) and how can they be obtained? 
To answer this question properly, it is necessary to examine 
first the nature of the Laplace variable s. This variable 
represents an operation in the time domain. S is a differentia- 
tion and 1/s is an integration. Thus it is intuitively felt 
that whatever f(z) is used, it should perform an approximate 
operation in the discrete time domain similar to that of the 
Laplace operator in the continuous time domain. Partarrieu 
[9] has shown that the selection of a particular function of 
Z to substitute for 1/s corresponds to the adoption of a 
discrete time numerical integration algorithm. A brief review 
of his development will be helpful. 


Consider a continuous transfer function 


3 Os gee a Ey (2-5) 


This transfer function represents the ratio of the transform 


ee the output, Y(s), to the transform of the input signal, X(s) 


H(s) = (2267) 





Be: 





This equation can be expressed as follows 











a ea) C=) ee (2-7) 
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Rearranging Eqs. (2-8) and (2-9) gives 


feel ; ec: 
b. s” W(s) = X(s) - £ b. s* W(s) (2=on 
: L 
1=0 ; 
m m 
Y(s) = £ a. s? Ws) (Gari 


Equations (2-10) and (2-11) can be depicted in an analog 
computer simulation as shown in Fig. 2.1. It can be seen that 
the replacement of the aru factors (integrators) by a function 
Geez denoted by [£(z)] 7 corresponds to a direct algebraic 
substitution for the variable s. If the hears is selected 
from the many numerical integration algorithms available in 
the literature [10], it follows that the substitution will, in 
fact, correspond to the operation of integration. How good 


a job it does, in the sense as to how well the discrete time 
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approximation represents the continuous seranstserenune-none 
depends upon which integration algorithm was selected. A 
study of several possible integration schemes and the effect 
of their approximations in the sinusoidal steady-state fre- 


quency domain forms the basis of this thesis. 


Be FREQUENCY RESPONSE CONSIDERATIONS 

A problem that the filter designer must face when formulat- 
mag a filter transfer function by the algebraic substitution 
method is that of translation from the S and the Z2 domain to 
the sinusoidal steady-state frequency domains. A common way 
to obtain the frequency response of a continuous filter is to 


let 
Ss = WwW (2=12)) 


and evaluate the magnitude squared and phase characteristics 
of the transfer function. For example, Suppose the desired 


transfer function is as shown 


H(s) = (2-13) 





where A(s) and B(s) are rational polynomials in s, then the 


magnitude function would be 


iz WV, 
pet = JAG) |" _ Aw) (2-14) 


IB(jw)]7 BB (w*) 


This expression can then be plotted and the filter classified 
as one of several types (low pass, band pass, band stop, etc.) 
An example of this procedure is shown in Fig. 2.2. 
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Transfer Function 





H(s) = = : = a<b 
; jw + a 
RA a 
S = jw 
: 2a ie + at ¥ eS 
(ju) |? = St 8, = Alas) 
W b Bl(w ) 
2 
|H (jo) | 
ab —_—_ a a ——s =e oP amr eel eee =~ 
2 
a 
be? HIGHPASS 
FOR a < b 
0 W 


Fig. 2.2. EXAMPLE OF FORMULATING FREQUENCY 
CHARACTERISTICS 2ROM CONTINUOUS 
TRANSFER FUNCTION. 
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A similar concept of filter frequency response exists for 
the digital filter. A designer will want to be able to examine 
the response from the specific filter and deternines ee 
frequency domain, it meets the needs for which it has been 


designed. To do this it is necessary to recall that [2] 


z= exp (sT) so that z= eJ (2-15) 


where T is the sampling period and s = jf. The symbol 2 is used 
here rather than w so as to distinguish the sinusoidal steady- 
state frequency for the continuous transfer function from the 
Sinusoidal steady-state frequency for the discrete time trans- 
fer function. Thus, to facilitate the formulation of the 
frequency response characteristics of the discrete transfer 
function directly from the sinusoidal steady-state character- 
istics of the continuous transfer function it iS convenient 


to have a second substitution 
WwW = FG) (2-16) 


that can be inserted in Eq. (2-14) to arrive directly at the 
magnitude squared frequency domain transfer function for the 
discrete case. A block diagram of this procedure is shown in 
Fig. 2.3. The derivation of these frequency-related functions 
erong with the resulting distortion of the filter character— 
metics of the transfer function is discussed in detail later 
Om. 

In the process of translating frequencies from one domain 


to another the designer encounters the question of whether 


ae 
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these frequency translations will be linear. In genenal they 
will not be linear and therefore a certain amount ofeatstogeren 
or "warping" of the frequency scales will result. This 
phenomenon poses certain restrictions and limitations on the 


use of the algebraic substitution as will be discussed. 


feo DLRECT PROCEDURE FOR FILTER SYNTHESIS 

The object of this section is to provide a detailed outline 
of the steps to be followed in employing the algebraic substi- 
futcion method for the direct synthesis of a digital filter: 

1. Determine the specifications. 

Before any actual design work can be done the designer 
must thoroughly understand his particular problem and ascertain 
what functions his filter must fulfill. The actual specifica- 
tions should be carefully studied and tabulated. These speci- 
fications may include cut-off frequencies, center frequencies, 
roll-off per octave, etc. The designer must have a clear 
picture of what his filter must do. 

2. Solve the approximation problem. 

This step encompasses the selection of the most suitable 
continuous filter that will approximate the desired filter. 
There are a wide variety of filters such as Butterworth, 
Chebyshev, Elliptic, and others that are described in detail 
ena their characteristics tabulated in various filter hand- 
books. The designer must be able to choose one of these 
classical filter designs to suit his particular problem and 


consequently he must be familiar with their characteristics. 
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3. Determine the transfer function. 

The transfer function of the selected continuous filter 
must be expressed as a rational function of S in sopders-omaeons, 
the algebraic substitution method. In addition the transfer 
function should also be expressed as a function of w (or ne) 
in order to have an insight into the frequency response of the 
fPieieer . 

4. Choose the numerical integration algorithm for the 

replacement of s. 

This step is probably the most important one in the chain. 
Whether the filter performs as it should depends to a great 
extent on the proper selection of the numerical integration 
algorithm. A familiarity with the various algorithms available 
for this purpose is essential. 
meeeeecrform the required algebra. 

Once the required substitution functions are known it is a 
relatively easy matter to accomplish the algebra required to 
mie the new G(z) in the form of the ratio of two polynomials 
in inverse powers of z. The designer may have to make 
allowances for the frequency distortion resulting from the 
non-linear frequency scale translation. The required frequency 
minmetiLon G(jf) can also be formed. This function will in 
general be fairly complex and will require effort for inter- 
pretation. 

6. Implement the transfer function. 
Once the necessary algebra has been accomplished and the 


transfer function has been formed it must be synthesized. 


AA. 





There are many references concerning the procedures to be 
employed in actually synthesizing tie sf ile map 

Be ean onan of this process are not within the purview 
of this thesis. Hess [11] has enumerated the twenty-four 
canonical forms for the synthesis of second-order sections 
and these may be either cascaded or paralleled to achieve the 
necessary implementation. 

7. Test the completed filter design. 

One method of observing the performance of the designed 
filter is to simulate the G(z) on a digital computer. Programs 
are available for this purpose. 

It can be seen from the foregoing that the critical points 
in the design procedure are incorporated in steps 4 and 5. 
The remainder of this thesis is directed to the study of 
integration algorithms and their application to the design - 


process. 
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Til. NUMERICAL ENTEGRA EC a 
SPECIFIC ALGORITHMS 

Numerical integration is the process of computing the value 
of a definite integral from a set of numerical values of the 
integrand. This process is accomplished by representing the 
integrand by some interpolation formula and then integrating 
this formula between the desired limits. The integrand is 
usually represented by a polynomial which approximates the 
actual shape of the function over the desired interval. Nearly 
all of the numerical integration methods employ this approx- 
imating technique and the basic difference between them is the 
complexity of the approximations. While, in general, the more 
complex an algorithm is, the more accurate it is; this does 
Hot mean that it is more desirable to use. This greater 
complexity also means more arithmetic operations to be performed 
and an increase in computation time. These may make a particu- 
lar algorithm unsuitable. 

This chapter is a study of three specific algorithms -- 
Euler methods, forward and backward, and the Trapezoidal 
method. All three can be expressed by the following discrete 


integration formula [12] 


Tk 

y* | (Ce) Gre (3-1) 
aaeea | 

Bu uke a kG apes GL sem.) BD) 
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where Y, and ¥,_, are the computed values of the integral at 
time uu and time ee respectively; f(T, ) is the value of the 


integrand at time T X is the parameter that determines the 


K* 
particular integration scheme. Consider Eq. (3-2) with the 


values of A = 0, l, = respectively. When 


A = 0 7 Vy = Yen] + T f(T.) Backward Euler (3-3) 

REN G MS Y.-7 + T f(T, _,) Forward Euler (3-4) 

Nes ive = ee Le Cee ec 
ye teil: ap k eal P 


These are the three algorithms and it can be seen that they 
@ieeter Only in how the input £(t) is handled. Fig. 3.1 asa 
representation of how these three algorithms approximate a 
function for integration. Each of the three methods will be 


discussed in the ensuing sections. 


A. FORWARD EULER INTEGRATION METHOD 

The Euler methods (both forward and backward) are also 
known as the rectangular rule. This stems from the geometric 
interpretation of the integration scheme. The integral of, or 
the area under, a curve can be approximated by rectangles as 
shown in Fig. 3.1. The area is brought up-to-date by adding 
the rectangle indicated to the area already accumulated. 


mos from Eq. (3-4) 


nr + T £(T ) (3-6) 


k kel jeeo Tt 


where f(T) _4) is the value of the function being integrated at 


time T The backward Euler as discussed in the next section 


Koa, * 
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.-l. GENERAL INTEGRATION FORMULA FOR 


SPECIFIC ALGORITHMS. 
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differs from the forward by using the present value of the 
aliqyonby ec f(T) to compute the present value of the output Yy- 
It can be seen that using the Euler method corresponds to 
approximating the input function by one that is piecewise 
constant. As a result the method is not very accurate when 
the function is changing very rapidly. In addition the 
integration time step will have to be small for the approxima- 
tion to have any validity. 

To derive the algebraic substitution for the forward Euler 


integration scheme, it 1s necessary to represent the discrete 


difference equation (3-6) in the z-domain. Recalling that 


z + = exp (-sT) eri) 


and that : 


Z [f(n-a)T]) = z°* 2Z[f(nT)] (3-8) 


memes possible to z-transform Eq. (3-6) yielding 


1 


Y(z) =z Y¥(z) + z+ 


te) (35) 


The integration transfer function is then 


(3-10) 





This transfer function representing integration in the s- 


domain is 


NG a = 
Pr (s) = (3-11) 


oR 
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A comparison of Eqs. (3-10) and (3-11) indicates the following 


correspondence can be made. 





(3-1) 


This then becomes the algebraic substitution for s for the 
forward Euler integration algorithm. 

It is now possible to examine the sinusoidal steady-state 
frequency domain characteristics of this s-plane to z-plane 
transformation. Letting s = jQ in Eq. (3-7) and applying this 


Pemag. (3-12) leads to 


= F3QT 
F(jn) = +S (3213) 
Te 
QT 
e eee Ee 
= £1 J sin >= (Baa 
F(jQ2) = a [sin QT + j(l-cos NT) } (3-15) 


This then is the required F(j2) for the forward Euler. This 
can be compared with the ideal integrator in the continuous 
sense to discover the frequency distortion which results. The 


magnitude and phase of Eq. (3-15) are 


2 sin un (3-16) 


[F(j2)| = 7 


ARG F(32) + 5 (3-17) 


These expressions can be compared to the magnitude and phase 


of the Laplace variable s when s = jw. 


a7 





> 7 


H(s) = s (3-18) 


H(jw) = jw (3-19) 
JH(jw) | = wo (2-310) 
ARG H(juw) = > | (3-21) 


Eqs. (3-16) and (3-20) can be combined to provide an analysis 


of the magnitude distortion between the continuous and discrete 


operations. 
TY Lee mee 9 0 6 _ 
ial eS a (3 22) 


This deviation from linearity for the magnitude as well as the 
phase angle deviation is shown graphically in Fig. 3.2. 
Although Eq. (3-22) is nonlinear it can be seen from Fig. 3.2a 
that there is a region of linearity close to the origin. MThus 
for small values of radian frequency the forward Euler magni- 
tude spectrum will closely approximate that of the ideal in- 
tegrator. This can be made more quantitative by considering 


the series expansion of Eq. (3-16) 





il 
FI} bh 
a, 

~ 

E 

~ 

WJ 
KH 
es 
++ 
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| F (52) | (3223) 


Therefore the magnitude spectrum of the forward Euler will 
closely resemble that of the ideal integration for values of 
2 such that 


Q3p2 


Q >> 74 


(3-25) 
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ideal integrator 





a) magnitude distortion 


ARG F(j2) 
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b) Phase angle distortion 


Fig. 3.2. FORWARD EULER FREQUENCY 
DISTORTION. 
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Ge equivalently 


¥24 


Q << a (3-26) 
V¥24 cies 
Q << 4 WS a ae : (3-27) 


The meaning of Eq. (3-27) can be made more explicit. Suppose 
there exists a continuous frequency spectrum whose highest 
frequency component is w+ In order for the Euler integration 
scheme to transform this spectrum into the discrete frequency 
domain without distortion, the sampling frequency (wo) must 
be much larger than 4/3 of the corresponding highest digital 
frequency component (2, ). If this condition is satisfied then 
the transformation will be linear. This implies that if a 
series of integrators are used in the continuous representation, 
then each transformation of these integrations to the digital 


domain must satisfy the criterion, namely, 


jae = © (3-28) 


Since T = 271/01 Eq. (3-28) puts a constraint on the value of 
the integration time step if a linear transformation is desired. 
However, there is a much more fundamental restriction which 
must be placed on T as the following discussion will show. 

The use of the Euler transformation should be limited to 
those continuous filters whose frequency spectrum is essentially 
band-limited. This is a very important result. Observe Fig. 


Ses which portrays the effect of the Euler integration scheme 
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Fig. 3.3. EULER TRANSFORMATION OF CONTINUOUS 
PREQUENCY SPECERUM) TO Discrete 
SebclLRUM: 
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on transforming a continuous filter spectrum into the digital 
filter spectrum. It can be seen that unless the continuous 
spectrum is band-limited the transformation will cause the 
high frequency ends of the spectrum to be "chopped-off" when 
being transformed to the digital representation. Since the 
transforming sine wave has finite amplitude the continuous 
Signal must have finite bandwidth. In addition the amplitude 
of the sine wave is directly related to this finite bandwidth 
in a manner analogous to the Sampling Theorem [13]. 

This theorem states that for a signal of finite bandwidth, 
Wh radians per second, a sampling rate equal to at least twice 
the maximum frequency (W) of the signal is necessary to recover 
the signal without distortion from an ideal low-pass filter. 
Fig. 3.4 represents this process. Fig. 3.4a shows a band- 


limited signal of W, radians per second. If the signal is: 


h 
sampled at a rate higher than that specified by the theorem, 
Fig. 3.4b is the result. However, if the signal is sampled 
at a rate less than the Nyquist rate, "aliasing" or "folding" 
@eeurs and this is illustrated in Fig.93.4c. 

To relate this to the Euler method, consider again Fig. 3.3. 
If the continuous filter has as its highest frequency component, 


Ws then to ensure that no part of this spectrum will be lost 


during the transformation, the value of T must be such that 


(e28)) 


As long as this condition is satisfied the amplitude of the 


transforming sine wave will be sufficiently large to guarantee 
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that all of the continuous spectrum will be mapped into the 
digital domain. This then is the fundamental constraint that 
LS Silvera on the value of T. This must always be satisfied 
when using the Euler method. However, it is necessary to set 
the value of T lower than that implied by Eq. (3-29) in order 
to keep the distortion reduced. 

Since the Euler method must be restricted to those contin- 
uous filters which have a band-limited spectrum, it is obvious 


that this technique is limited. 


B. BACKWARD EULER INTEGRATION METHOD 

The backward Euler, as was pointed out in the previous 
section, differs from the forward Euler algorithm in the way 
the input is treated. The backward method uses the present 
value of the input to compute the present value of the output. 


Re@all Eq. (3-3) 


Y, = ¥,_, + T £(7,) (3229) 


and it can be seen that the input and output are considered at 
the same time point. Intuitively this would appear to be more 
accurate and would not involve any prediction on the part of 
the integration scheme, and in general this is true. 

An alternative way to understand the difference between the 


two Euler methods is to look at the differential equation 


oe Ee) (2=an 
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The Euler method integrates this type of vditierent al sean aeen 


by approximating the derivative as shown 
y — ea oo (3=32) 


In the forward method the derivative 1S approximated at time 
T._, and the formula of Eq. (3-4) resulles but ehes haem ais 


Euler approximates the derivative at time T, and Eq. (3-3) is 


k 
applicable. When selecting an algorithm it must be decided 
whether or not the present value of the input will be available 
to compute the present value of the output. It would be pre- 
ferable to use the present value of the input, but there can 
be circumstances which make this impossible. 

Similar to the procedure followed in the previous section 
the z-transform can be ssoiied to Eq. (3-30) "and the sresule 


is 


Pi) = bs Sele) So sale Bane 


From this equation the integration transfer function can be 


femme cd . 





= (2) = (2-34) 


Comparing this to Eq. (3-11) yields 


owe 
Ss +> = (3-35) 





for the backward Euler integration algorithm. To find the 


Sinusoidal steady-state response of Eq. (3-35) the 


S15. 





substitutions mentioned previously are applied. This results 


in the following 


1 4OT 
F(j2) = == (Gaea 


F(j2) = 2 [sin of - j(1-cos 97)] Con 


Comparing this expression to Eq. (3-15), it is seen that only 
the sign of the real part has changed. This does not affect 
the magnitude of F(j®2) which is the same as Eq. (3-16). Thus 
all of the analysis that was presented concerning the frequency 
distortion for the forward Euler is true for the backward 
Euler and the same constraints on the value of the integration 
time step are applicable. 

The phase angle characteristics of the backward case, 
however, are different from the forward Euler and are given 


by 


ARG F(j2) = > - (3-38) 


ae 
b 


This is plotted in Fig. 3.5 along with the magnitude distor- 
tion. It can be seen that the phase angle of the backward 
Euler iS approximately the same as the ideal case when 


= >> — Thus for low-pass filters with low cutoff frequencies 
the phase angle of the backward Euler is virtually the same as 
the ideal integration. 

One property that would be desirable for any integration 


scheme to possess would be the mapping of the imaginary axis 
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b) Phase angle distortion. 


Fig. 3.5. BACKWARD EULER FREQUENCY 
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of the s-plane on to the unit circle of )the Z-plancis tiem mes 
circle in the z-plane establishes the region of Stabile, seen 
the Z-domain; namely, the interior of this unit circle. If an 
algorithm mapped the imaginary axis of the s-plane onto the 
unit circle of the Z-plane, then it would be possible to 
ensure that as long as the continuous filter was stable then 
the digital representation would also be stable. 
To examine this characteristic closer consider the back- 
word Euler integration transfer function Eq. (3-35). “Muilltzpiie 


numerator and denominator by z yields 








1 el ie 
Ss = T z (3 39) 
Solving for z 
dl z 
4 * [-st (3-40) 
Evaluating this expression with s = jw 
_ _1 = 


it is obvious that the magnitude of z does not equal unity. 


ime fact 
1 
|z| = ——>—- (3-42) 
Gagne \ Oe 
ARG z = fan ~ WT (3-43) 


The above mapping pair transforms the imaginary axis of the 


S-plane to a pair of circles that are completely enclosed by 
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the unit circle. This is plotted in Figg 6b.) Tomemenounse 
that the backward Euler will map continuous system poles 

enees in the left half of the s-plane into the stable region 
of the z-plane. The forward Euler, however, does not perform 
in this manner as shown in Fig. 3.6a. The appropriate equations 


for the forward scheme are 


lz) = 1 + w*2?)l/ (3-44) 


AReee = tone. our (a=e5) 


It is seen that the imaginary axis is mapped onto a curve which 
goes to infinity along asymptotes of + 90 degrees. Thus with 
the forward Euler algorithm there iS no assurance that a pole 
in the left half of the s-plane will be mapped into the stable 
region. | 

The foregoing analysis would seem to indicate that the 
backward Euler integration would be superior to the forward 
algorithm for use in the algebraic substitution method. 
Although they possess the same frequency magnitude characteris-— 
tics the phase angle characteristics of the backward Euler are 
Superior to that of the forward. Certainly as far as stability 
is concerned the backward scheme is superior. The next sec- 


tion discusses the trapezoidal integration algorithm. 


C. TRAPEZOIDAL INTEGRATION METHOD 
The trapezoidal method of discrete integration derives its 


name from the geometric shape of the approximation that is 
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made to the function being integrated. Instead of using a 
piecewise constant approximation (rectangles) as was done by 
the two Euler methods, the trapezoidal scheme uses a trapezoid 
to fit the function as shown in Fig. 3.7. This results in a 
more accurate fit to the actual function and reduces the error 


between the actual and computed values. Recalling Eq. (3-5) 
Y, = ¥,_, + = [£(t,) + £(T,_,)] (3.46) 
k k-1 Z k k-1 ° 


it can be seen that for the trapezoid algorithm the solution 


to the differential equation 


y = f(t) (2a4a) 


1s being approximated at the average of the two time points 
Ty and Toye This characteristic ensures a "good fit" for’ 
functions which are montonically increasing or decreasing. 
Also this algorithm is quite good if the input function is 
reasonably well-behaved. 

To obtain the algebraic substitution for the trapezoidal 


method it is necessary to z-transform Eq. (3-46) and doing 


this yields 


(Se = we > (F(z) Be lai (3.48) 


The integration transfer function can now be formed and is 





Mi _ T ltz Za 
7 (z) = 5 (3-49) 
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Fig. 3.7. TRAPEZOIDAL APPROXIMATION TO 
AREA UNDER CONTINUOUS CURVE. 
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An examination of Eq. (3-49) reveals that the integration 
transfer Function possesses a bilinear nature. In fact several 
authors [l, 5] refer to the trapezoidal transformation as the 
bilinear substitution. This bidimearte aa. one property of 
this algorithm which makes it very useful for digital filter 
design. Proceeding as before, it is possible to make the 


correspondence between the s-plane and the z-plane 


ae ca coos (3-50) 


This is the algebraic substitution for the trapezoidal inte- 
gration scheme. 

By algebraic manipulation of Eq. (3-50), one ot the mesr 
desirable features of the bilinear transformation can be 
demonstrated. It was stated earlier that it would be convenient 
for a transformation to map the imaginary axis of the s-plane 
onto the unit circle of the z-plane. The trapezoidal algorithm 
possesses this property. To see this consider Eq. (3-50) and 


disregard the constant factor 2/T. 


—~ 274 E 
Ser ior] Ca 


Selving for z yields 





ye l+s (3250) 
l-s 
Let s = jw 
2 1+j¥ 7 jz] = 1 (3-53) 
eget 
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This transformation is the only one which possesses this 
property of mapping the imaginary axis onto the unit circle 
without increasing the order of the system. This ensures that 
a pole in the left half of the s-plane in the continuous 
system will be mapped into the stable region of the z-plane. 
This fact makes the trapezoidal integration algorithm a very 
useful tool for the algebraic substitution method. 

To examine the frequency domain characteristics of the 
trapezoidal method, apply the substitutions mentioned in 


Section A. to Eq. (3-50). The result is 


=92T 
2 1-e J 
F(j2) = 7 Tae aT (3-54) 
F(j@) = 3% tan = R55) 


The first item of interest to be noted is that the F(j2) has 
no real part and, consequently, the phase of the trapezoidal 
algorithm is precisely the same as that of the ideal operator. 
This characteristic is very important when developing the 
frequency substitutions as will be seen in Chapter IV. To 


see the magnitude distortion compare Eqs. (3-19) and (3-55) 


jw = j é tan (5) 
= = ban _ Geom 


This magnitude distortion is clearly shown in Fig. 3.8. It 


can be seen from Fig. 3.8 that in the low frequency case the 
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Fig. 3.8. TRAPEZOIDAL ALGORITHM FREQUENCY 
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transformation function can be considered linear and, in fact, 


his transformation is linear for values of 2 such that 


1 
Oa 5 w (3-58) 


The above condition is similar to that given for the Euler 
transformations in the previous section. Both conditions result 
from the fact that the trigonometric terms involved are linear 
about the origin. However there is an important difference 
between the trapezoidal and the Euler methods which is also 
attributable to these trigonometric terms. Fig. 3.9 shows 
the transformation from a continuous spectrum to a digital 
spectrum. The fact that the trapezoidal transformation 
contains the tangent factor means that the entire continuous 
frequency band can be transformed. In fact, the entire 
Semernuous Spectrum 1s transformed Inte the antervalli ae 
in the discrete domain. This property makes the bilinear 
substitution a more universal type of transformation as it 
can be used for any spectrum. 

Unfortunately, for the filter designer, there still exists 
the problem of distortion when using the trapezoidal method. 
This distortion can be overcome, in part, by pre-warping [5] or 
frequency scaling the particular continuous realization so that 
after the algebraic substitution is made the selected critical 
frequencies are at the desired values. However this will work 
only for those frequencies selected and may yield an unsatis- 
factory design at other frequencies. Another solution and one 


that has been mentioned previously would be to increase the 
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sampling frequency, thereby ensuring that the critical 

frequencies would be much less than the sampling frequency. 
This solution also has drawbacks. Increasing the sampling 
frequencies also reduces the time available for computation. 
If the designer is faced with a fixed computing speed then 
increasing the sampling frequency would limit the complexity 
of the filter than can be implemented. Thus reducing the 
sampling interval might not be a desirable tool to use. In 
this case pre-warping might be the only feasible solution. 

mane trapezoidal or bilinear transtermaeton sas Peenmiced 
more extensively than any other substitution technique due to 
its universal application to the frequency domain and the fact 
that it does not increase the complexity of the filter; l.e., 
an 2am order continuous filter is transformed into an ee order 
discrete filter. Also the fact that the trapezoidal algorithm 
possesses ideal phase angle characteristics greatly simpiifies 


the frequency substitution as will be demonstrated in Chapter 


IV. 
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IV. NUMERICAL INTEGRATION - 
GENERAL ALGORITHM 
Chapter III presented the general case discrete integration 


formula and it is repeated here for emphasis 
eae) Vy act mk f(T, _4) ge ie Er a=) 


The two quantities in this formula that are of particular 
interest are T (the integration time step) and A. The parameter 
A may be considered a weighting parameter in that the selection 
of A’ determines which time point of f(t) will have the greatest 
effect on the integration process. These two quantities (T 

and A) determine the type of performance that can be achieved 
uSing this integration formula. It was seen earlier that 

= and 1 led to three well- 


known integration algorithms. However there is no absolute 


choosing the values of i to be 0, 


restriction on this value of X and there may be instances 
where another choice would provide better results. 


It is known, however, that a value of A such that 
O< A <5 (4-2) 


guarantees stability of the numerical integration process. 

This constraint on the value of } means that as long as the 
poles of the continuous system lie in the left half of the 

s-plane, the solution to the system will be numerically 


stable for any positive value of T. If X > = then the value 
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of T will have some upper bound that will depend on the exact 
value of A that is used and on the location Gt the wroctouormede 
Spans Appendix B contains a detailed analysis of these 
stability conditions. It is sufficient hereto benravareuc= 

the advantages that accrue from choosing A eeoecans to 

Eq. (4~2). 

Fig. 4.1 is a representation of the general integration 
formula for various values of X. It can be seen that the 
approximation to the actual function moves up or down as A 
is varied. Thus as A is made smaller the approximating 
rectangle grows larger and, conversely, as dX is made larger 
the approximating rectangle becomes smaller. Thus it is seen 
‘that the value of 4» determines which time point will have the 
greatest effect on the approximation. 

It would be convenient to have the general integration 
formula available for use in the algebraic substitution method. 
This would greatly facilitate implementing the method once the 
designer has selected the parameters he wants to use. This 
chapter contains the derivation of the general case transfor- 
mations that could then be employed once the values of T and 
X have been selected. In addition the derivation of the 
general case frequency substitutions are presented in section 
Be. 

A. GENERAL ALGORITHM TRANSFORMATIONS 
Following the procedures outlined previously it 1s 


necessary to z-transform Eq. (4-1) to put it in the form 
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required for the algebraic substitution. Doing this yields 


y(z) = 27> y(z) + ta 2) F(z) + T(1-A) F(Z) (4-3) 


Solving to form the integration transfer function produces 


= IL 


Y(Z) (1-2 =) = TIA@ 7a (4-4) 

Y TL (Z-~-1) A421] 

FF (2) = = (4-5) 
C7) 


Eq. (4-5) is the desired integration transfer function and 


from this the following correspondence can be made. 


“1 
4) (4-6) 


ee: 
Mit, ST ay 


This then is the required algebraic substitution for s for the 
general integration algorithm. By selecting a value for i 

one can immediately obtain the desired substitution for s. 

7 1, Eq. (4-6) becomes in turn 
moee( 3-35), Eq. (3-50), and Eg. (3-12). 


It can be seen that when A = O, 


The next transformation required is that for mapping from 
the continuous frequency domain to the discrete frequency 


domain. To accomplish this consider Eq. (4-6) and let 


i = exp (-j2T) 


Applying this to Eq. (4-6) 


jQT 
=jNT 


leu 


[1+A (e 


F(jq) = # (4-8) 


=) 


a2 





# 


This form can be expressed in terms of trigonometric quantities 


as follows 
sin OT 
2 


Pie@) = o2 (4-9) 


OT : elk 
cos 57— -j C2A—1) es aa > 


This form permits the necessary F(j2T) expression to be quickly 
developed for substitution into the frequency domain transfor- 
mations which are derived in the next section. 

To analyze these generalized results it is helpful to 
rearrange Eq. (4-9) into a form similar to the ones used in 


the preceding chapter. Eq. (4-9) can be written 


> J sin ut cos = ee IL) sin? 
es oT 2 OF eo ame 
cos” 5— +) (2 X= ese = 
a... : QT 
Dividing through by sin 5 yields 
Oye = 4 eeu 22 
a ; 2 2 a (4-11) 
(2A-1) <2 (elone ee 
Consider the magnitude of both sides of Eq. (4-11). 
2 1 
LT at. eye (4-12) 
T feot? S24 (24-1) °1777 


NO 


Eq. (4-12) permits an examination of the frequency distortion 
that will be present when this general integration algorithm 
is used. 

Fig. 4.2 is a representation of Eq. (4-12) for various 


values of A. It can be seen that when A = 0, , and 1 the 
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plots are the same as in the previous chapter. However, for 
other values of X} there are no clear-cut analytical expres- 
sions for the curves. The amplitude of the curves increase 

aS X approaches = So that the trapezoidal rule (A = =) is the 
limiting upper bound. Conversely, the amplitude of the curves 
decrease as A approaches either 0 or 1 so that the Euler 
methods are the limiting lower bound. 

In the previous chapter it was discovered that the Euler 
transformations should be used only for band-limited spectra. 
Also the value of T that was used would have to be constrained 
in such a manner as to ensure that all of the spectrum would 
be transformed without loss. Here in the general case it can 
be seen that the same situation exists. All of the transfor- 
mations except for the trapezoidal one have finite amplitudes 
and thus should only be used for finite spectra. HEC, 
even though the amplitudes are finite they can be made arbi- 
trarily large by the proper choice of X. This means that the 
constraint on the value of T may be relaxed. Unfortunately, 
linearity considerations must also be taken into account and 
these will cause the value of T to be kept small to avoid 
excessive distortion. 

Another factor that must be considered at this point is 
just how complex a transformation one would want to use. It 
4S obvious from Eq. (4-6) that a value of A other than 0, * 
and 1 will make the algebra associated with that particular 
substitution more cumbersome. In addition any value of xX 
other than 1 will cause some phase distortion and this may be 


Z 
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an important consideration. In fact, section B will show that 
when the phase is not ideal, the frequency transformations 


become much more involved. 


B. GENERAL CASE FREQUENCY TRANSFORMATION 
It was pointed out in Chapter II that it would be convenient 


to have a substitution 
W = F(IL) (4-13) 


that could be used to directly obtain the magnitude squared 
representation for the discrete transfer function. For each 
integration algorithm that has been considered a F(j2) has 
been derived. Now if this expression was the proper substitu- 
tion for w then this problem would not only be solved, but 
greatly simplified as these expressions were not very eho 
However, in general, these expressions do not fulfill the 
task. To clearly demonstrate this, a simple example will be 
presented. 


Consider the transfer function depicted in Fig. 2.2 


E 
Olle (s) _ s + a (4-14) 





If the backward Euler transformation is applied to this 
continuous expression, a transfer function in the z-domain 


momopbtained. 





E male 
out a atatiz— (aE) 
in l+btT+z 





edit and obtain the magnitude squared function. 
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2 a? 
— (42) » (l + aT cos — +. aan OT (4-16) 
in (1 + b T =) cos Gmie sine 


This is the correct expression as it was derived directly from 
the discrete transfer function. However this procedure would 
usually be very time-consuming and a better solution would be 
to directly substitute an expression into the magnitude squared 


Sinusoidal steady-state transfer function. 


S 
cout (5 
un 


we 
= ay (4-17) 








However when this is done using the F(j2) derived earlier the 


Pome tLon obtained is 
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It is evident after only a little algebra that Eq. (4-16) and 
Eq. (4-18) are not equivalent, and consequently the procedure 
employed is not valid in this case. 

The negative result of the above example demonstrates the 
need for a general transformation which will be valid for 
every case. Two separate cases were considered in developing 
these general transformations. The first case was the linear 
term (s + a) aS was used in the preceding example. The second 


was the quadratic factor (s? + as + b). Since transfer function 


a) 4) 





factors should be put into no more than second order terms 
for ease in implementing the filter it was considered sufficient 
to derive general transformations for these two cases cnly. 


Consider the linear term situation. 





(s + a) = [£(z) + all, Ae (4-19) 
Now examine the magnitude squared function. 
ne +h ac = ees a5 alle ao) = F(el@Ty (4-20) 
= (a + Pe fel eee et (4-21) 
Deeee i= a> Gagne F(e 2") + Rae. oii Coy 
ue = DE We (ei, raed) (i (4-23) 


mimes then is the final result for the Pinear case. Tt Ganvbe 


seen that the term (2a Re ples is the one which prevents 


the simpler transformation 


2 


2 
Wo = IP(es Fy | 


(4-24) 
from being valid in the general case. For those transforma- 
tions which do not have a real part, e.g., trapezoidal rule, 
the general transformation of Eq. (4-23) reduces to the much 
eeneber form of Eq.em(4—-24) 9 It can alse be seen taat wien a 
transformation does. have a real part, this is equivalent to 
Saying that the phase angle characteristics differ from those 
of the ideal integrator. This is the case with both Euler 


methods. 
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The case of the quadratic factor requires considerable 
mathematical manipulation and yields a transformation which is 


quite a bit more complex. Assuming a quadratic factor of the 


pmobaltt| 


2S 57 +as tb (4-25) 


the magnitude squared function becomes 


Blac) = weet uc (ac - OR (4-26) 


By a procedure Similar to that followed for the linear case 


the following transformation was obtained 


w4 ts w* (a*-2b) + b* = Bee) | a 7° a’ Face joes 
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ja Rel o(el eo ve 
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) + 2b R, Fe (eI Ty 


+ 2a Deke Eee (aon 


Thus the entire factor in the continuous domain becomes that 
shown in Eq. (4-27) in the discrete domain. Even in this case, 
however, for those algorithms which possess ideal phase 
characteristics the general transformation reduces to a much 
Simpler and more logical form. Unfortunately, in general, a 
Simple substitution will not suffice. 

These transformations provide a direct avenue to forming 
the frequency domain representations of a digital filter. 
While they do not appear simple in form, in many cases they will 


reduce to a very convenient expression which will provide the 
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designer with the desired information regarding the frequency 
response of a particular filter. They will also provide a 
good deal of insight into the frequency distortion that will 


be present. This is demonstrated in the next section. 


eee FE REOURNCY DISTORTION CONSIDERATIONS 

It has been shown that every algorithm that can be derived 
from the general integration formula possesses some degree of 
meecuency distortion. In Chapter III this distortion was 
explicitly analyzed for three specific algorithms. In addi- 
tion guidelines were presented for minimizing this distortion. 
These rules indicated that by decreasing the sampling interval 
(T) appropriately, the amount of distortion can be reduced to 
a negligible value over a specified frequency range. In 
general this is the case with any algorithm that is used. As 
long as the sampling frequency is much greater than the cut- 
off frequency the transformation involved will be essentially 
linear. However, the general integration formula also involves 
the parameter A and certainly it would be helpful to know which 
value of X is optimum for reducing distortion. 


To examine this consideration further recall Eq. (4-12) 


il 


Z 
QoQ = — (4-28) 
T foot* S24 (ja-1) 71°"? 
This equation can also be written as 
2 ean So 


fl + (2A-1) caine ea10/? 
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At this point it is necessary to apply series expansion 
approximations. These approximations are valid as the function 
here is being examined in the vicinity of the origin and the 
argument (5) 1s very small. Hence higher order terms can be 


mrscounted. This leads to 


3 
io het 
w 2 2 — (4-30) 
[1 + (20-1)? (5% 4 GEN 4) 1/2 
w = 2 OE y OM. 2 (ayeay2 Be D742) 4-31) 
w LY 2 2 2 6 
: 1 ~ a , 
Since 172 et ) provided syst IL 
La) 
rien 
22 or 91 |.) 2 eee : 
w = 2 SEs (-b (aaa)? + 4H SY) (4-32) 


At this point it is possible to form the ratio of the nonlinear 
term and the linear term to obtain a representation of the 


fractional deviation from linearity. 


(a) = (4233) 


NE 


Fractional Deviation = A = [5 - 


The item of interest that can be seen from Eq. (4-33) is that 
the value of A will be small since it involves a second order 
term. Also the factor [= = 5 (20-17) ] can be analyzed for 

different values of A and the result is very interesting. The 


value of A that gives the largest first order distortion 


oF 





1s = and there are two values that result in no first order 
distortion. These are 0.09 and 0.91. ObvicuSsly cnemcmlaenee 
be interested in the larger value since it would not insure a 
stable integration. However, the smaller value would insure 
this stability while not contributing any first order 
aestortion. 

While this analysis is only valid in the vicinity of the 
Origin it must be remembered that the designer will want to 
adjust the sampling interval so that the approximation will be 
in this linear zone. The important point of this analysis is 
that the designer will be faced with decisions about which 
algorithm to use and it may be necessary to make allowances 
for a certain amount of distortion in order to Nave a less 
complex transformation. However, the analysis that has been 
presented will permit the designer to make intelligent cHeaieee 
concerning the necessary parameters in employing the algorithms. 
At this point it would be instructive to look at an example 
of a filter design problem to gain insight into how the theory 
applies in practice. 

Example. The continuous transfer function mentioned in the 
previous section can be used to understand the method of 
attacking the approximation problem. The magnitude squared 


function of that filter was 


2 wy" + a’ 
H(w') = ———s (a < bd) (4-34) 
2 2 
Ww + b 
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Also the discrete counterpart can be formed using the transfor- 


mation derived in section B., 


2 02 
2 - 
Wein (lot oa T = 6S eect mL (4-35) 


(l + b T -—- cos Qn) ¢ + Stun” Cor 


This can be written as follows. 


=. (l + a T) Sune se a’ 
, 2 Aly 

7) (i+ b LT) sean Sn b 

At 


Fig. 4.3a shows a plot of both Eq. (4.34) and Eq. (4-36) where 
Eq. (4-36) is merely shown in general form for any arbitrary 
value of T. 

Since this particular continuous spectrum is not band- 
limited, the Backward Euler transformation was not the speared 
one to use. However, its use illustrates the ideas that are 
involved and also show that part of the spectrum is lost as 
a result of the transformation. Fig. 4.3a clearly shows the 
periodicity of the filter frequency response. In addition it 
is interesting to note that if T is allowed to approach zero, 
the digital filter response approaches that of the continuous 
filter. This is of course consistent with the theory of 
sampling. This also confirms the assertion made earlier that 
for very small values of the sampling interval the transforma- 
tion is essentially linear. 

Thus it is possible to reproduce the continuous filter 


frequency response as close as necessary if the value of T 
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is adjusted properly. For this example if one wished to have 
the digital filter response approximate that of the continuous 
filter for values of w up to ten then a value of T such that 
T = 0.01 will do the job. This is shown in a magnified view 
in Fig. 4.3b. Similarly any such specification can be 
Satisfied if the designer is willing to make the value of T 
small enough. 

The transformations covered in this chapter offer the 
designer a variety of approaches to the design problem. While 
one transformation may achieve greater linearity than another, 
the second may be easier to implement. Therefore, the filter 
designer is faced with the recurrent engineering problem. 

What selection of the available choices can be put together 
to give the best possible design? The analyses presented will 


make the decisions less difficult. 
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V. SUMMARY AND SUGGESTIONS 
FOR FURTHER RESEARCH 

The important points presented in this thesis are: 
1) The algebraic substitution method is a convenient and 
relatively simple procedure to use in designing a digital f2keer 
from continuous filter characteristics. While the discrete 
filter is only an approximation to the in eee one, it is 
possible to make this discrete realization approach the 
desired realization as close as necessary. 
2) The problem of frequency distortion is present ineany OL 
the transformations that would be used in the algebraic 
substitution technique. However, the effect of this distoreLlem 
can be offset by intelligent application of the guidelines 
contained in this thesis. These guidelines include the rules 
for achieving the desired amount of linearityodn. thie = trans Forma. 
tions, the stability constraints, and the general case frequency transformations. 

The foregoing studies indicate that an interesting area 
for further exploration and analysis would be a computer 
study of these various substitution techniques employing the 
theory presented herein. This would involve the use of 


different algorithms with different classical filters. 
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APPENDIX A 


STABILITY CONSIDERATIONS 


An important aspect of any numerical integration method 
is the stability of the method [14]. In this context Sa ersvepil dL 1 ey 
refers to whether or not the errors generated at each step one 
the numerical integration process tend to decay (stable) or 
grow (unstable) as the solution progresses in time. The 
question of stability becomes important in the selection of how 
small an integration time step must be used. It would be very 
convenient for an integration algorithm to be stable for all 
positive values of T as long as the eigenvalues of the system 
being processed are in the left half of the s-plane. 

To develop a general stability criterion consider the 


general algebraic substitution. 





dl 
S = = 1? (Agee) 
(z “-1)A+l 
Solving for Zz. 
_ Seabee 

2 = STA-STTL (A.2) 
Now substitute a root at Ss = O+tju. 

je = ATo + ATjo + 1 (on) 


~ TotjwT(A-l) + 1 
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Bor 0Stabidi ty (2c sor 





Z 25202 
ot IA : 
eS oul (A. 4) 
(1+T95)\-To)~ + (wAT-wT) 
Breer much algebra the result is: 
226 
TE LEESIPA 2) cee 7a (Az) 
OD th 


For a root in the left half of the s-plane the value of o will 
Besnegative. Therefore the quantity on the right becomes a 
negative number. The quantity on the left will be greater 


than or equal to zero as long as 
0O< AS (A.6) 


Therefore as long as Eq.(A.6) is satisfied then any positive 
value of T will ensure the inequality in Eq.(A.5) is satisfied 
and the integration will be stable. This does not imply that 
the solution will be unstable for i > > , but rather that a 
constraint on the size of T must be imposed to ensure stability. 
This is precisely the case with the forward Euler. However, 


the Trapezoidal and backward Euler are completely stable in 


mrs Sense. 
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